package scu.maqiang.homogenization;

import static scu.maqiang.numeric.Constants.GP2;
import static scu.maqiang.numeric.Constants.GW2;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3H201;
import scu.maqiang.fes.FES3H81;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.Pair;
import scu.maqiang.numeric.SRMatrix;

public class SOTS3DH8201 {
	public SOTS3DH8201(FES3H81 fsH8, FES3H201 fsH20) {
		this.fsH8 = fsH8;
		this.fsH20 = fsH20;
	}

	/**
	 * ��������ȴ��������ͷ����ȴ���ϵ����һ�ס����׵������Լ����Ȼ��ȴ���ϵ����
	 * 
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 * @param SOCS
	 */
	public void KxxCellSolution(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D,
			double[][] SOCS) {
		SRMatrix A = new SRMatrix(fsH8.GetNdof());
		fsH8.assembleStiff(kxxFunc, param, BVPType.COMMON, A);
		double[][] RHS = new double[3][fsH8.GetNdof()];
		fsH8.assembleGlobalVector(kxxFunc, param, this::elementkxxFirstOrder, BVPType.COMMON, RHS);
		fsH8.applyBC_MBN(A, Direct.All, 1, 2, 3, 4, 5, 6);
		fsH8.applyBC_MBN(RHS[0], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
		fsH8.applyBC_MBN(RHS[1], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
		fsH8.applyBC_MBN(RHS[2], Direct.All, 0.0, 1, 2, 3, 4, 5, 6);
		NewIterSSolver solver = new NewIterSSolver(A);
		MVO.fill(FOCS, 0.0);
		solver.PCGSSOR(RHS[0], FOCS[0], 1.5, 1);
		solver.PCGSSOR(RHS[1], FOCS[1], 1.5, 1);
		solver.PCGSSOR(RHS[2], FOCS[2], 1.5, 1);

		computeHomogenizedKxx(kxxFunc, param, FOCS, D);

		RHS = new double[9][fsH8.GetNdof()];
		fsH8.assembleGlobalVector(kxxFunc, param, FOCS, D, this::elementKxxSecondOrder, BVPType.COMMON, RHS);
		MVO.fill(SOCS, 0.0);
		for (int i = 0; i < 9; i++) {
			fsH8.applyBC_MBN(RHS[i], Direct.X, 0.0, 1, 2, 3, 4, 5, 6);
			solver.PCGSSOR(RHS[i], SOCS[i], 1.5, 1);
		}
	}

	/**
	 * �����ȴ���ϵ��һ�׵��������Ҷ���������
	 * 
	 * @param xyz
	 * @param coef
	 * @param tp
	 * @param eleXYZ
	 */
	public void elementkxxFirstOrder(int[] ele, double[][] xyz, double[][] coef, BVPType tp, double[][] eleXYZ) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] v = new double[8];
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		Arrays.fill(eleXYZ[0], 0.0);
		Arrays.fill(eleXYZ[1], 0.0);
		Arrays.fill(eleXYZ[2], 0.0);
		for (int i = 0; i < 2; i++) {
			double xi = GP2[i];
			for (int j = 0; j < 2; j++) {
				double et = GP2[j];
				for (int k = 0; k < 2; k++) {
					double zt = GP2[k];
					double detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
					double c = MVO.dot_product(coef[0], v);
					double wxarea = GW2[i] * GW2[j] * GW2[k] * detJ * c;
					for (int m = 0; m < 8; m++) {
						eleXYZ[0][m] += -dx[m] * wxarea;
						eleXYZ[1][m] += -dy[m] * wxarea;
						eleXYZ[2][m] += -dz[m] * wxarea;
					}
				}

			}
		}
	}

	/**
	 * ��һ�׵������������ȴ������Ȼ�ϵ������
	 * 
	 * @param kxxFunc
	 * @param param
	 * @param FOCS
	 * @param D
	 */
	public void computeHomogenizedKxx(ScalarFunc[] kxxFunc, double[] param, double[][] FOCS, double[][] D) {
		double[][] xyz = new double[8][3];
		Mesh3H8 mesh = (Mesh3H8) fsH8.getMesh();
		double[] NN1 = new double[8];
		double[] NN2 = new double[8];
		double[] NN3 = new double[8];
		double wholevolume = 0.0;
		double[][] coef = new double[kxxFunc.length][8];
		MVO.fill(D, 0.0);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		double[] v = new double[8];
		double xi, et, zt;
		for (int idx = 0; idx < mesh.getNt(); idx++) {
			int[] ele = mesh.getElement(idx);
			mesh.getCoordPerElement(ele, xyz);
			MVO.getSlice(FOCS[0], ele, NN1);
			MVO.getSlice(FOCS[1], ele, NN2);
			MVO.getSlice(FOCS[2], ele, NN3);
			double[] x = MVO.col(xyz, 0);
			double[] y = MVO.col(xyz, 1);
			double[] z = MVO.col(xyz, 2);
			fsH8.getCoefFromFunc(kxxFunc, xyz, mesh.getElementLabel(idx), param, coef);
			for (int i = 0; i < 2; i++) {
				xi = GP2[i];
				for (int j = 0; j < 2; j++) {
					et = GP2[j];
					for (int k = 0; k < 2; k++) {
						zt = GP2[k];
						double detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						double mKxx = MVO.dot_product(coef[0], v);
						double wxarea = GW2[i] * GW2[j] * GW2[k] * detJ * mKxx;
						D[0][0] += (1 + MVO.dot_product(NN1, dx)) * wxarea;
						D[0][1] += MVO.dot_product(NN2, dx) * wxarea;
						D[0][2] += MVO.dot_product(NN3, dx) * wxarea;
						D[1][0] += MVO.dot_product(NN1, dy) * wxarea;
						D[1][1] += (1 + MVO.dot_product(NN2, dy)) * wxarea;
						D[1][2] += MVO.dot_product(NN3, dy) * wxarea;
						D[2][0] += MVO.dot_product(NN1, dz) * wxarea;
						D[2][1] += MVO.dot_product(NN2, dz) * wxarea;
						D[2][2] += (1 + MVO.dot_product(NN3, dz)) * wxarea;
					}

				}
			}
			wholevolume += CG.volume3DH8(xyz);
		}
		MVO.mul(D, 1.0 / wholevolume);
	}

	/**
	 * �����ȴ���ϵ�����׵����������Ҷ���������
	 * 
	 * @param xyz
	 * @param kxx
	 * @param FOCS
	 * @param homo
	 * @param tp
	 * @param eleXYZ
	 */
	public void elementKxxSecondOrder(int[] ele, double[][] xyz, double[][] kxx, double[][] FOCS, double[][] homo, BVPType tp,
			double[][] eleXYZ) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] v = new double[8];
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		for (int i = 0; i < 9; i++) {
			Arrays.fill(eleXYZ[i], 0.0);
		}
		for (int i = 0; i < 2; i++) {
			double xi = GP2[i];
			for (int j = 0; j < 2; j++) {
				double et = GP2[j];
				for (int k = 0; k < 2; k++) {
					double zt = GP2[k];
					double detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
					double mkxx = MVO.dot_product(kxx[0], v);
					double wxarea = GW2[i] * GW2[j] * GW2[k] * detJ;
					for (int kk = 0; kk < 8; kk++) {
						eleXYZ[0][kk] += ((mkxx - homo[0][0] + mkxx * MVO.dot_product(dx, FOCS[0])) * v[kk]
								- mkxx * MVO.dot_product(v, FOCS[0]) * dx[kk]) * wxarea; // 11
						eleXYZ[1][kk] += ((-homo[0][1] + mkxx * MVO.dot_product(dy, FOCS[0])) * v[kk]
								- mkxx * MVO.dot_product(v, FOCS[0]) * dy[kk]) * wxarea; // 12
						eleXYZ[2][kk] += ((-homo[0][2] + mkxx * MVO.dot_product(dz, FOCS[0])) * v[kk]
								- mkxx * MVO.dot_product(v, FOCS[0]) * dz[kk]) * wxarea; // 13
						eleXYZ[3][kk] += ((-homo[1][0] + mkxx * MVO.dot_product(dx, FOCS[1])) * v[kk]
								- mkxx * MVO.dot_product(v, FOCS[1]) * dx[kk]) * wxarea; // 21
						eleXYZ[4][kk] += ((mkxx - homo[1][1] + mkxx * MVO.dot_product(dy, FOCS[1])) * v[kk]
								- mkxx * MVO.dot_product(v, FOCS[1]) * dy[kk]) * wxarea; // 22
						eleXYZ[5][kk] += ((-homo[1][2] + mkxx * MVO.dot_product(dz, FOCS[1])) * v[kk]
								- mkxx * MVO.dot_product(v, FOCS[1]) * dz[kk]) * wxarea; // 23
						eleXYZ[6][kk] += ((-homo[2][0] + mkxx * MVO.dot_product(dx, FOCS[2])) * v[kk]
								- mkxx * MVO.dot_product(v, FOCS[2]) * dx[kk]) * wxarea; // 31
						eleXYZ[7][kk] += ((-homo[2][1] + mkxx * MVO.dot_product(dy, FOCS[2])) * v[kk]
								- mkxx * MVO.dot_product(v, FOCS[2]) * dy[kk]) * wxarea; // 32
						eleXYZ[8][kk] += ((mkxx - homo[2][2] + mkxx * MVO.dot_product(dz, FOCS[2])) * v[kk]
								- mkxx * MVO.dot_product(v, FOCS[2]) * dz[kk]) * wxarea; // 33
					}
				}

			}
		}
	}

	public void assembleHeatSOTSSolution(FES3H81 fsFine, double[] T_0, double[] T_1, double[] T_2, FES3H201 fsHomo,
                                         double[] T_Homo, FES3H81 fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		double[] xyz = new double[3];
		double[] cellXYZ = new double[3];

		double[][] grad = new double[3][fsHomo.GetNdof()];
		double[][] Se = new double[6][fsHomo.GetNdof()];
		// System.out.println("Computing the gradient of the homogenized solution!...");
		// fsHomo.computeGradient2(T_Homo, grad);
		// fsHomo.computeSecondDerivative(T_Homo, Se);
		System.out.println("Computing the first and second derivative of the homogenzied solution!...");
		fsHomo.computeFirstAndSecondDerivative(T_Homo, grad, Se);
		double val = 0;
		double corrector1, corrector2;
		System.out.println("Assembling started!...");

		double[][] homoCoord = new double[20][3];
		double[][] homoCoef = new double[20][20];
		int[] ele = fsHomo.mesh.getElement(0);
		int[] eleIdx20 = ele.clone();
		double[] eleVal20 = new double[20];
		fsHomo.mesh.getCoordPerElement(ele, homoCoord);
		CG.H20ShapeFunctionCoef(homoCoord, homoCoef);
		Pair<Double, Integer> homoData = new Pair<>(0.0, 0);

		double[][] cellCoord = new double[8][3];
		double[][] cellCoef = new double[8][8];
		ele = fsCell.mesh.getElement(0);
		int[] eleIdx8 = ele.clone();
		double[] eleVal8 = new double[8];
		fsCell.mesh.getCoordPerElement(ele, cellCoord);
		CG.H8ShapeFunctionCoef(cellCoord, cellCoef);
		Pair<Double, Integer> cellData = new Pair<>(0.0, 0);

		for (int i = 0, n = fsFine.GetNdof(); i < n; i++) {
			fsFine.mesh.getCoordPerNode(i, xyz);
			cellXYZ[0] = xyz[0] * N - (int) (xyz[0] * N);
			cellXYZ[1] = xyz[1] * N - (int) (xyz[1] * N);
			cellXYZ[2] = xyz[2] * N - (int) (xyz[2] * N);

			fsHomo.getValue(T_Homo, xyz, homoData, homoCoef, eleIdx20, eleVal20);
			T_0[i] = homoData.getFirst();

			fsCell.getValue(FOCS[0], cellXYZ, cellData, cellCoef, eleIdx8, eleVal8);
			val = cellData.getFirst();

			corrector1 = 0.0;
			corrector1 += val * MVO.dot_product(grad[0], eleIdx20, eleVal20);
			corrector1 += MVO.dot_product(FOCS[1], eleIdx8, eleVal8) * MVO.dot_product(grad[1], eleIdx20, eleVal20);
			corrector1 += MVO.dot_product(FOCS[2], eleIdx8, eleVal8) * MVO.dot_product(grad[2], eleIdx20, eleVal20);
			T_1[i] = T_0[i] + ep * corrector1;
			corrector2 = 0.0;
			corrector2 += MVO.dot_product(SOCS[0], eleIdx8, eleVal8) * MVO.dot_product(Se[0], eleIdx20, eleVal20); // xx
			corrector2 += MVO.dot_product(SOCS[1], eleIdx8, eleVal8) * MVO.dot_product(Se[3], eleIdx20, eleVal20); // xy
			corrector2 += MVO.dot_product(SOCS[2], eleIdx8, eleVal8) * MVO.dot_product(Se[5], eleIdx20, eleVal20); // xz
			corrector2 += MVO.dot_product(SOCS[3], eleIdx8, eleVal8) * MVO.dot_product(Se[3], eleIdx20, eleVal20); // xy
			corrector2 += MVO.dot_product(SOCS[4], eleIdx8, eleVal8) * MVO.dot_product(Se[1], eleIdx20, eleVal20); // yy
			corrector2 += MVO.dot_product(SOCS[5], eleIdx8, eleVal8) * MVO.dot_product(Se[4], eleIdx20, eleVal20); // yz
			corrector2 += MVO.dot_product(SOCS[6], eleIdx8, eleVal8) * MVO.dot_product(Se[5], eleIdx20, eleVal20); // xz
			corrector2 += MVO.dot_product(SOCS[7], eleIdx8, eleVal8) * MVO.dot_product(Se[4], eleIdx20, eleVal20); // yz
			corrector2 += MVO.dot_product(SOCS[8], eleIdx8, eleVal8) * MVO.dot_product(Se[2], eleIdx20, eleVal20);// zz
			T_2[i] = T_1[i] + ep * ep * corrector2;
			System.out.println(i + "   of  " + n + " has been assembled!");
		}
	}

	public void assembleHeatSOTSSolution_Cube(FES3H81 fsFine, double[] T_0, double[] T_1, double[] T_2,
                                              FES3H201 fsHomo, double[] T_Homo, FES3H81 fsCell, double[][] FOCS, double[][] SOCS, int N, double ep) {
		double[] xyz = new double[3];
		double[] cellXYZ = new double[3];

		double[][] grad = new double[3][fsHomo.GetNdof()];
		double[][] Se = new double[6][fsHomo.GetNdof()];
		// System.out.println("Computing the gradient of the homogenized solution!...");
		// fsHomo.computeGradient2(T_Homo, grad);
		// fsHomo.computeSecondDerivative(T_Homo, Se);
		System.out.println("Computing the first and second derivative of the homogenzied solution!...");
		fsHomo.computeFirstAndSecondDerivative(T_Homo, grad, Se);
		double val = 0;
		double corrector1, corrector2;
		System.out.println("Assembling started!...");

		double[][] homoCoord = new double[20][3];
		double[][] homoCoef = new double[20][20];
		int[] ele = fsHomo.mesh.getElement(0);
		int[] eleIdx20 = ele.clone();
		double[] eleVal20 = new double[20];
		fsHomo.mesh.getCoordPerElement(ele, homoCoord);
		CG.H20ShapeFunctionCoef(homoCoord, homoCoef);
		Pair<Double, Integer> homoData = new Pair<>(0.0, 0);

		double[][] cellCoord = new double[8][3];
		double[][] cellCoef = new double[8][8];
		ele = fsCell.mesh.getElement(0);
		int[] eleIdx8 = ele.clone();
		double[] eleVal8 = new double[8];
		fsCell.mesh.getCoordPerElement(ele, cellCoord);
		CG.H8ShapeFunctionCoef(cellCoord, cellCoef);
		Pair<Double, Integer> cellData = new Pair<>(0.0, 0);

		for (int i = 0, n = fsFine.GetNdof(); i < n; i++) {
			fsFine.mesh.getCoordPerNode(i, xyz);
			cellXYZ[0] = xyz[0] * N - (int) (xyz[0] * N);
			cellXYZ[1] = xyz[1] * N - (int) (xyz[1] * N);
			cellXYZ[2] = xyz[2] * N - (int) (xyz[2] * N);

			fsHomo.getValue_Cube(T_Homo, xyz, homoData, homoCoef, eleIdx20, eleVal20);
			T_0[i] = homoData.getFirst();

			fsCell.getValue_Cube(FOCS[0], cellXYZ, cellData, cellCoef, eleIdx8, eleVal8);
			val = cellData.getFirst();

			corrector1 = 0.0;
			corrector1 += val * MVO.dot_product(grad[0], eleIdx20, eleVal20);
			corrector1 += MVO.dot_product(FOCS[1], eleIdx8, eleVal8) * MVO.dot_product(grad[1], eleIdx20, eleVal20);
			corrector1 += MVO.dot_product(FOCS[2], eleIdx8, eleVal8) * MVO.dot_product(grad[2], eleIdx20, eleVal20);
			T_1[i] = T_0[i] + ep * corrector1;
			corrector2 = 0.0;
			corrector2 += MVO.dot_product(SOCS[0], eleIdx8, eleVal8) * MVO.dot_product(Se[0], eleIdx20, eleVal20); // xx
			corrector2 += MVO.dot_product(SOCS[1], eleIdx8, eleVal8) * MVO.dot_product(Se[3], eleIdx20, eleVal20); // xy
			corrector2 += MVO.dot_product(SOCS[2], eleIdx8, eleVal8) * MVO.dot_product(Se[5], eleIdx20, eleVal20); // xz
			corrector2 += MVO.dot_product(SOCS[3], eleIdx8, eleVal8) * MVO.dot_product(Se[3], eleIdx20, eleVal20); // xy
			corrector2 += MVO.dot_product(SOCS[4], eleIdx8, eleVal8) * MVO.dot_product(Se[1], eleIdx20, eleVal20); // yy
			corrector2 += MVO.dot_product(SOCS[5], eleIdx8, eleVal8) * MVO.dot_product(Se[4], eleIdx20, eleVal20); // yz
			corrector2 += MVO.dot_product(SOCS[6], eleIdx8, eleVal8) * MVO.dot_product(Se[5], eleIdx20, eleVal20); // xz
			corrector2 += MVO.dot_product(SOCS[7], eleIdx8, eleVal8) * MVO.dot_product(Se[4], eleIdx20, eleVal20); // yz
			corrector2 += MVO.dot_product(SOCS[8], eleIdx8, eleVal8) * MVO.dot_product(Se[2], eleIdx20, eleVal20);// zz
			T_2[i] = T_1[i] + ep * ep * corrector2;
			System.out.println(i + "   of  " + n + " has been assembled!");
		}
	}

	public void assembleHeatSOTSSolution(FES3H81 fsFine, double[] fineSol, FES3H201 fsHomo, double[] T_Homo,
                                         FES3H81 fsCell, double[][] FOCS, double[][] SOCS, int N, double ep, double[][] xyzs, double[][] value) {

		double[] eleVal20 = new double[20];
		double[] eleFineVal8 = new double[8];
		double[] eleCellVal8 = new double[8];
		double[][] grad = new double[3][fsHomo.GetNdof()];
		double[][] Se = new double[6][fsHomo.GetNdof()];
		System.out.println("Computing the first and second derivative of the homogenzied solution!...");
		// System.out.println("Computing the gradient of the homogenized solution!...");
		// fsHomo.computeGradient2(T_Homo, grad);
		fsHomo.computeFirstAndSecondDerivative(T_Homo, grad, Se);

		double val = 0;
		double corrector1, corrector2;
		double[] cellXYZ = new double[3];

		// int previousFineEleIdx = 0;
		double[][] fineCoord = new double[8][3];
		double[][] fineCoef = new double[8][8];
		int[] ele = fsFine.mesh.getElement(0);
		int[] eleFineIdx8 = ele.clone();
		fsFine.mesh.getCoordPerElement(ele, fineCoord);
		CG.H8ShapeFunctionCoef(fineCoord, fineCoef);
		Pair<Double, Integer> fineData = new Pair<>(0.0, 0);

		// int previousHomoEleIdx = 0;
		double[][] homoCoord = new double[20][3];
		double[][] homoCoef = new double[20][20];
		ele = fsHomo.mesh.getElement(0);
		int[] eleHomoIdx20 = ele.clone();
		fsHomo.mesh.getCoordPerElement(ele, homoCoord);
		CG.H20ShapeFunctionCoef(homoCoord, homoCoef);
		Pair<Double, Integer> homoData = new Pair<>(0.0, 0);

		// int previousCellEleIdx = 0;
		double[][] cellCoord = new double[8][3];
		double[][] cellCoef = new double[8][8];
		ele = fsCell.mesh.getElement(0);
		int[] eleCellIdx8 = ele.clone();
		fsCell.mesh.getCoordPerElement(ele, cellCoord);
		CG.H8ShapeFunctionCoef(cellCoord, cellCoef);
		Pair<Double, Integer> cellData = new Pair<>(0.0, 0);

		for (int i = 0, n = xyzs.length; i < n; i++) {
			double[] xyz = xyzs[i];
			cellXYZ[0] = xyz[0] * N - (int) (xyz[0] * N);
			cellXYZ[1] = xyz[1] * N - (int) (xyz[1] * N);
			cellXYZ[2] = xyz[2] * N - (int) (xyz[2] * N);
			// val = fsFine.getValue(fineSol, xyz, eleIdx8, eleVal8);
			fsFine.getValue_Cube(fineSol, xyz, fineData, fineCoef, eleFineIdx8, eleFineVal8);
			value[i][0] = fineData.getFirst();
			// previousFineEleIdx = fineData.getSecond();

			// val = fsHomo.getValue(T_Homo, xyz, eleIdx20, eleVal20);
			fsHomo.getValue_Cube(T_Homo, xyz, homoData, homoCoef, eleHomoIdx20, eleVal20);
			value[i][1] = homoData.getFirst();
			// previousHomoEleIdx = homoData.getSecond();

			// val = fsCell.getValue(FOCS[0], cellXYZ, eleIdx8, eleVal8);
			fsCell.getValue_Cube(FOCS[0], cellXYZ, cellData, cellCoef, eleCellIdx8, eleCellVal8);
			val = cellData.getFirst();
			// previousCellEleIdx = cellData.getSecond();

			corrector1 = 0.0;
			corrector1 += val * MVO.dot_product(grad[0], eleHomoIdx20, eleVal20);
			corrector1 += MVO.dot_product(FOCS[1], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(grad[1], eleHomoIdx20, eleVal20);
			corrector1 += MVO.dot_product(FOCS[2], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(grad[2], eleHomoIdx20, eleVal20);

			value[i][2] = value[i][1] + ep * corrector1;
			corrector2 = 0.0;
			corrector2 += MVO.dot_product(SOCS[0], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(Se[0], eleHomoIdx20, eleVal20); // xx
			corrector2 += MVO.dot_product(SOCS[1], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(Se[3], eleHomoIdx20, eleVal20); // xy
			corrector2 += MVO.dot_product(SOCS[2], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(Se[5], eleHomoIdx20, eleVal20); // xz
			corrector2 += MVO.dot_product(SOCS[3], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(Se[3], eleHomoIdx20, eleVal20); // xy
			corrector2 += MVO.dot_product(SOCS[4], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(Se[1], eleHomoIdx20, eleVal20); // yy
			corrector2 += MVO.dot_product(SOCS[5], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(Se[4], eleHomoIdx20, eleVal20); // yz
			corrector2 += MVO.dot_product(SOCS[6], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(Se[5], eleHomoIdx20, eleVal20); // xz
			corrector2 += MVO.dot_product(SOCS[7], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(Se[4], eleHomoIdx20, eleVal20); // yz
			corrector2 += MVO.dot_product(SOCS[8], eleCellIdx8, eleCellVal8)
					* MVO.dot_product(Se[2], eleHomoIdx20, eleVal20); // zz
			value[i][3] = value[i][2] + ep * ep * corrector2;
			// System.out.println("solutions in " + Arrays.toString(xyz) + ": " +
			// Arrays.toString(value[i]));
			System.out.println(i + "   of  " + n + " has been assembled!");
		}
	}

	FES3H81 fsH8;
	FES3H201 fsH20;
}
